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Abstract 

We have made substantial advances in elucidating the properties of the susceptibil- 
ity of the square lattice Ising model. We discuss its analyticity properties, certain closed 
form expressions for subsets of the coefficients, and give an algorithm of complexity 
0(A^^) to determine its first N coefficients. As a result, we have generated and analyzed 
series with more than 300 terms in both the high- and low-temperature regime. We 
quantify the effect of irrelevant variables to the scaling-amplitude functions. In partic- 
ular, we find and quantify the breakdown of simple scaling, in the absence of irrelevant 
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scaling fields, arising first at order \T — Tc\^f^ , though high-low temperature symmetry 
is still preserved. At terms of order |r — Td^''/^ and beyond, this symmetry is no longer 
present. The short-distance terms are shown to have the form (T — rc)^(log \T — Tc\Y 
with p > q^. Conjectured exact expressions for some correlation functions and series 
coefficients in terms of elliptic theta functions also foreshadow future developments. 

Keywords 

Ising susceptibility, high-temperature series, low-temperature series, scaling function, irrele- 
vant variables, differentiably finite functions, scaling fields. 

1 Introduction 

Since Onsager's [|^ celebrated solution of the Ising model free energy in 1944, followed by 
Yang's proof of Onsager's result for the spontaneous magnetization in 1952, almost half 
a century has passed during which time many, if not most of the world's most able mathe- 
matical physicists have devoted themselves to the problem of elucidating the susceptibility. 
Attempting to list all these contributions would produce a bibliography of prohibitive length, 
and one that would inevitably commit many sins of omission. Therefore rather than attempt 
this, we will only make mention of those papers that have directly motivated our work here, 
and crave the forgiveness of those who we have inadvertently offended. 

While much of the notation for describing the square lattice Ising model is standard, 
we begin by defining our notation here both for the benefit of the more casual reader and 
to emphasize those cases where we deviate from convention. The interactions in the two 
perpendicular directions are taken to be 

K = I3J, K' = (3J'. (1.1) 

but we also often set K' = K to discuss the isotropic lattice. For high temperatures, 
s = sinh2i^ and s' = sinh2i^' are appropriate variables for series expansions 0, while for 
low temperatures, we use 1/s and 1/s' instead. Thus, 

s* = sinh2K* = l/sinh2K = 1/s, s'* = sinh2Js:'* = l/sinh2K' = 1/s'. (1.2) 

In many cases, high-temperature and low-temperature formulas can be obtained from each 
other by Kramers- Wannier duality with a simple interchange of primes and stars. The critical 
temperature is defined by the condition s' = s*. 
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A conventional high-temperature variable is f = tanh K, while an often-used low-temper- 
ature variable is m = exp(— 4/^). The translations between these and our variables are 

s = smk2K = 2v / {I - v"^) (1.3) 

and 

s* = 2u^/^l{l-u), (1.4) 

and similarly for the primed variables. 

In studying the critical behavior, we will use both the variable t = 1 — T^/T and more 
frequently 

r = (1/s — s)/2 (isotropic) (1.5) 

to parameterize deviations from the critical temperature. To leading order, r = 2X^^/21. 
An elliptic parameterization will be useful, and to that end we define the elliptic modulus, 

{s'/s* = ss' for T>T„ 
k=l ' (1.6) 

ys*/s' = l/ss' for T<Tc. 
Let (Tjj- be the spin at lattice site (^, j) and define the two-point function 

C(M,iV) = (ao,oaM,iv). (1.7) 
In terms of the correlation functions the susceptibility is 

where M. is the magnetization. 

In 1956 Syozi and Naya |^ presented an approximation to the anisotropic high-tempera- 
ture susceptibility which gave the correct critical point, correct critical exponent, an am- 
plitude estimate that was wrong by less than 1%, and reproduced the first 8 terms of the 
series expansion. It was also exact along the disorder lineQ In 1976 a celebrated paper by 
Wu, McCoy, Tracy and Barouch showed how the high- and low-temperature expansions 

^For the fully anisotropic triangular lattice Ising model, with coupling constants Vi = tanh Ji/fcT, i — 
1,2,3, the condition for the disorder line is that viV2 + = 0. Along this line the correlations decay 
exponentially, and the partition function factorizes |^ . The disorder condition is of no particular interest in 
the case of the nearest neighbor square lattice Ising model. 
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of the susceptibihty could be understood in terms of a multi-particle expansion, with an 
odd number of particles being appropriate at high temperatures and an even number at low 
temperatures. With this interpretation it became clear that the result of [^] was just the 
lowest order, or one-particle approximation to the susceptibility. 

In terms of the elliptic modulus k ( |1.6| ) the high- and low-temperature susceptibilities 
can be written 

oo 

i3-\^ = k-Hi-k'y^J2x^''^'^ (1-9) 

and 

oo 

p-'x- = {l-k'f*Y.^^''^ (1.10) 
1=1 

respectively, where x*^"'^ is the sum over all lattice separations of the j-particle contribution 
to the two-point function, and was first given as a 2j-fold multiple integral. It was 
subsequently shown that this integral can be reduced to a j-fold integral of the form 

x(^) = y ■ ■ ■ y du,{G'^^^ff^\ (1.11) 

where G^^^ is a fermionic determinant and f^^^ is an algebraic function. This reduction 
has been achieved by various routes i,0-[ll. The factor, G^^\ appearing in the integrand, 
which can be expressed in terms of Pfaffians 0, has been found to have a product form, 
first by Palmer and Tracy in the low temperature regime [Q, and then independently by 
Yamada in both the high- and low-temperature regimes f^, |10|- From the product form 
it readily follows that the first non-zero term in (|LTT|) is 2^^^'i)y^/^. From the original 



expression one could only conclude that each integral entered at order fc-'/^, so clearly massive 
cancellations occur. While this comes as a surprise if handling the integral directly, it 
follows straightforwardly [Tl|] from the product form. 

In terms of the elliptic modulus, the first terms in the high- and low-temperature expan- 
sions for the isotropic {K = K') susceptibility are 



k 
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xLo = 7-—T7: (1-12) 

(1 — A;2 

and 

(,)_ (1 + P)E-(1-P)K 

37r(l-A:)(l-A;2) ' ^ ' 

where E and K are the complete elliptic integral^ of the second and first kind respectively. 
^We trust there is no confusion with the couphng constant K. 
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Anisotropic versions of these formulae are given in section ^. Simple forms for higher terms 
in the expansion are not known. 

In [|12| one of us gave compelling evidence that unlike the free-energy and spontaneous 
magnetization, the anisotropic susceptibility K') is not differentiably finite. A series 
in n variables, /(z), is said to be differentiably finite or D-finite if and only if it satisfies a 
system of n partial differential equations of the form 

d d^* 



P,,o(z)/(z) + P,i(z)— /(z) + ■ ■ ■ + P,,^(z)^/(z) = 0, (1.14) 



where the Pj,j(z) are polynomials and for each i = 1, ■ ■ ■ ,n, Pi^ki (z) is not the null polynomial, 
see e.g. Proposition 2.2 in |13]. Thus the expression for the susceptibility was shown to be 
in a different — and less tractable — class of function than other known properties of the Ising 
model. The evidence for this was based on the observation (not proved) that the anisotropic 
susceptibility x(^?'^^')? a function of v with v' fixed has a natural boundary on the unit 
circle |f | = 1. 

For the isotropic susceptibility, another of us 0,11] provided strong confirmation (though 
again, not a proof) of this observation by showing that the circle |s| = 1 in the complex 
s = sinh2ir plane is a natural boundary. These two observations are discussed further in 



section |37^ , following a discussion of the general anisotropic case in section [O 



In section p.3| we also prove the important result that while x(^) is not D-finite, x^-'\k) 
is D-finite for all j. 

Two other directions in which we have achieved substantial progress are in the generation 
of series coefficients for the individual series x^^\ continuing work initiated in [§, [ill, and 
even greater progress in obtaining the coefficients of the total series x using nonlinear partial 
difference equations for the correlation functions |ll^- |T^] . 

In order to generate the series for the total susceptibility x+ or X- without comput- 
ing separately the j-particle contributions, a more efficient method of series generation is 
obtained by first returning to the expression ( |1.8| ) of the susceptibility as the sum over all 
lattice separations of the two-point correlation functions. 

In the scaling limit, the two-point functions were found to satisfy a nonlinear differential 
equation of Painleve type which was then solved to give the leading scaling terms |t|~^/^ 
and ItI"^/"^ exactly [H,|l^. In 1980, a discrete analogue of this equation was discovered by 
McCoy and Wu fl^ which holds for the two-point functions of the lattice Ising model at 
arbitrary temperature, and which reduces to the Painleve equation in the scaling limit. In 



the same year, a simple set of partial difference equations was derived by Perk [|T^], which 
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reproduced the equation of McCoy and Wu and provided an additional equation. Also in 
the same year, an unrelated set of difference equations, which can be used to compute the 



correlation functions on the diagonal, M = N, was obtained by Jimbo and Miwa [jl6| . Many 
of these developments are described in some detail in 

The difference equations are valid for arbitrary temperature and were used by Kong, et 
al. to obtain exactly the leading "short distance" constant terms in the susceptibility 
both at the ferromagnetic and anti-ferromagnetic points at T = ±Tc. This work was later 
extended to give the amplitudes of the term r log |r| [^,|2l|]. Here we dramatically extend 
that work by obtaining all terms in the "short-distance" part| -Bf/af (see (|1.16|) and (|1.17|) ) 
of the susceptibility to 0(r^^). Details necessary for the generation of C{M,N) appear in 
section ^ while in section |6.1| we describe the numerical analysis of the C{M, N) that leads 
us to conclude that the "short-distance" terms have the form 

BiM = E E ^i/af ^'(log \r\r (isotropic). (1.15) 

g=0 p=0 

Although clearly nonanalytic at r = we have denoted these "short-distance" terms in 
( |1.15| ) by i?, as a further reminder that they also include the analytical "background." The 
actual coefficients in ( |1.15| ) can be found in the Appendix. 

The quadratic difference equations of Perk [|T^ can also be used to generate high- and 
low-temperature series for x a^nd as shown in section § the series coefficients can be obtained 
in polynomial time! 

While some people have expressed the view that a polynomial time algorithm for the 
computation of the series coefficients constitutes a solution, it is clearly preferable to have a 
closed form expression. Nevertheless, a polynomial time algorithm is equivalent to a complete 
solution if one seeks only the series coefficients, as to expand any closed form expression also 
takes polynomial time. As the history of the development of these key nonlinear recurrences 
described above shows, the ingredients for such an algorithm have existed unexploited in the 
literature for many years. 

Analysis of the resulting series of hitherto unimaginable length combined with the "short- 
distance" knowledge contained in ( |1.15D leads to a solidly based conjecture specifying com- 
pletely the remaining "scaling" part of the susceptibility x of the isotropic Ising model. Near 
the anti-ferromagnetic point the "short-distance" terms are complete and we simply have 
for T>Tr 



•^When we speak of the "short-distance" part, we mclude the analytic background term, often ignored in 
scahng discussions of the critical region. 
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P ^Xaf = 5af (isotropic). (1.16) 
Near the ferromagnetic point we conjecture (for T > Tc or T < Tc) 

f3~\± = Co^{2K,V2y/^\T\-'^/^F± + B{ (isotropic) (1.17) 

where the scahng-amphtude functions F± have (possibly asymptotic) expansions in integer 
powers of r without any of the powers of log |r| present in -Bf/af- The leading terms are 

F± = 1 + r/2 + 5rV8 + 3rVl6 - 23rV384 - 35rV768 + /jfV^ + O(r^) (isotropic) 

(1.18) 

where f^^ ^ In fact the breakdown in equality is dramatic, and we estimate that 
ff^ = -0.1329693327... , and ff^ = -6.330746944... , where more accurate values of 
these and further terms in (|1.18|) through order r^^ are given in the Appendix. In section |6.2| 
we describe the "short-distance" subtraction and analysis on which the assumed form of the 
expansion of the scaling-amplitude function F± (|1.18|) is based, while the analysis leading to 
the numerical values of the coefficients in ( |1.18| ) is sketched in section |6.3| . 



Aharony and Fisher |2^, ^ have predicted a scaling-amplitude function F(A&F) that 
is equal above and below on the assumption that the Ising model critical region can be 
described entirely by two nonlinear scaling fields. Our exact result ( p.. 181 ) is clearly different 
and furthermore the explicit expansion (cf. eqs. (22-24) in |Tl|) 



r') (1.19) 



F(A&F) = 1 + r/2 + 5rV8 + 3r7l6 - llrVl92 - 17x^384 + 97r73072 + 0( 

differs from ( |1.18|) at order r^. This is unequivocal evidence for the presence of at least one, 
and almost surely two, irrelevant operators^. There is further possible evidence for irrelevant 
operators in the "short-distance" terms ( p..l5| ) which contain powers of log|r| beyond the 
first starting at T^(log |t|)^, and thus are not of the "energy" form given by the nonlinear 
field analysis. 

An important numerical study investigating corrections to scaling was that of Gartenhaus 
and McCullough who confirmed the F(A&F) form in ( |1.19| ) through O(r^) and provided 



Note that r changes sign in ( 1.18| ) and ( |1.19| ) as we change from T > Tc to T < Tc- 
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a good estimate of the term linear in r in Bf in ( |1.17|) . Estimates of terms to 0(r log |r|) 



in i?af in ( |1.16| ) were obtained in |^3[. An attempt |Tl| to go beyond this using longer series 



than available in was inconclusive other than to indicate the necessity of terms beyond 



that predicted by Aharony and Fisher P^,pB 



Our numerical work began in part as a modest attempt to improve on |TT[ but expanded 



to where it now clearly quantifies the effect of nonlinear scaling fields as well as irrelevant 



operators. As anisotropy is a marginal operator, extending the present calculation to the 
anisotropic square lattice would, we expect, be extremely helpful in better understanding 
the effect of the irrelevant operator(s) we have identified. Both this and a study of the 
susceptibility on the triangular and hexagonal lattices are projects we hope to tackle in the 
near future. 

The layout of the paper is as follows. In section ^ we define the model and give some 
useful parameterizations. In section ^ we show how the key integrals referred to above may 
be simplified, and provide a short proof of the assertion that x^^\ so defined, is D-finite, 
even though, as we have seen, x is presumably not. In section ^ the computation of the 
susceptibility series from the correlation functions by means of nonlinear partial difference 
equations is shown to be achievable in polynomial time. In section |^ we first discuss the 
isotropic series in k, and then the q series form of the susceptibility. In this subsection 
some regularity features of the coefficients are discovered and the consequences partially 
developed. In section ^ we summarize our numerical work, state our conjecture giving the 
complete analytic structure of the isotropic susceptibility, and quantify the effect of irrelevant 
variables in the scaling fields. Finally, in section |^ we review scaling theory as it applies to 
the two-dimensional Ising model, and comment on the relevance of our results to this theory 
and to the renormalization group. The high- and low-temperature series for Xiso are given 
on the WWW at site www.ms.unimelb.edu.au/~tonyg. 

2 Definitions and notation 

In this section of the paper it will be convenient to formulate the general anisotropic model. 
Our later numerical work is confined to the isotropic model. 



The complementary modulus of the elliptic modulus k (p.. 61) is given by k' = yl^^W 



Here the prime is not related to the anisotropy. The moduli, k and fc', are related to the 
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elliptic nome, g, by 

The connection between the two moduli yields a well-known theta function identity. (See 
section 13.20 of Bateman for these and other formulas.) 

In terms of these variables the magnetization takes a particularly simple form P7[] 

oo _ 2n— 1 

M = {l-{ss'rY^ = l^\ = \{^^^ (2.2) 

n=l 

for T < Tc, and otherwise. 

The theta functions introduced above have useful infinite sum and product forms, which 
we will subsequently use. For zero argument these are 

oo oo 

^2(0, g) = 2 5^ = 2gi/^ J](l - q'-){l + q'^f, (2.3) 

71=0 n=l 
oo oo 

^3(0, g) = 1 + 2 5^ g"^ = - + (2.4) 

n=l n=l 

oo oo 

^4(0, g) = 1 + 2 $^(-1) = 11(1 - g2")(l - g2"-i)2. (2.5) 



n=l n=l 



Thus, 



1 I 2n '-^ 1 2 



n=l 



n=l 



In the elliptic parameterization, one variable is taken to be either k or g. We take the 
other variable to be the anisotropy parameter. Onsager |^ used the variable a defined by 

iis* forT>r„ 
snm = < (2.7) 
I is' for T < Tc. 

We will also need the related variable 

, ii/s' ioTT>T,, 
sma' = { ' (2.8) 
i/s* for T < Te. 
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In some contexts, it will be useful, following ref. [|T^, to use instead the variables a and a' 
defined by 



cot a 
cot a' 



— zV^snia 
—i\/ksnia' 



V s's* 



(2.9) 
(2.10) 



Obviously cot a = tana'. The isotropic values of these variables are snia = snm' = i/Vk 
and a = a' = n/A. 

We observe that g^/^ = x where x is the variable of ref. |27|, defined by 



-2K 



X 



8n-7 



)(1 



X 



n=l 



X 



8n-5 



)(1 



X' 



8n-3^ ' 



(2.11) 



in which the spontaneous magnetization has a simple product form. The elliptic modulus of 
ref. [27 1 is related to our A; by a Landen transformation [^, section 15.6]. 

In reference [^| it was also noted that the first terms of the isotropic {K = K') high- 
temperature susceptibility have a simple product expression, which breaks down at order 
g,8/4 This is explained by the expansion of the high-temperature susceptibility in multi- 
particle states described in the Introduction. The contribution of one-particle states has a 
product form, and the first three-particle state contributes at order g^/^. This one-particle 
product form is 



^(1) 

Also 



,1/4 -Q (i + g"/T(i-g^")^ 



,1/4 



ii (l-gn/4)2(l + ^n)2 ^3(0, g) 



n 

n=l 



l + g"/^)2(l- 



2nN3 



;i - g"/4)2 



(2.12) 



(Strictly speaking the product form above is not the product form of ref. which actually 
breaks down only at order g^/^. They differ by a factor which first contributes at order g^/^ 
and the absence of this factor in ref. exactly compensates for the addition of the first 
term of xLlM 

Attempts to do the same for the low-temperature series have failed as we do not know 

^ (2) 

of a product form for xiso- The function K has one, namely 

K=|7re2(0,g), 

while the function E is given in terms of g by the formula |2^ 

er(0,g) 



(2.13) 



^ g|(0,g) + gKO.g) i. 1 



3^1(0, g) 



12K ^;(0,g)- 



(2.14) 
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3 Integral formulae for x^^^ 

In this section we present integral expressions for the x^^^ and show that they define D-finite 
functions. 

3.1 Trigonometric form of integrals 

We will start with the form given by Yamada [P, [To|, p9H3T[| for the j-particle contribution 



X^^^ in terms of elliptic variables, and derive from it an expression in terms of trigonomet- 
ric/hyperbolic variables. The integral expression is 



where the modulus of the complete elliptic integral K is fc, and this same modulus is assumed 
in all Jacobi elliptic functions which appear below. Considered as a function of one of the 
Uni the integrand has a single simple pole on the real axis which derives from the last factor 
in the integrand. The contour of integration is deformed in the vicinity of the pole so that 
only half the residue is taken. The function G^^"^ can be written as the product 

G(^'^= n (3.2) 

l<m<n<j 

and z„, Xn and hmn are given by 

iu)n \ ~'~ ^'^'^ '^'^ 



sn I (un — ia') sn ia' dn u„ — dn ia' sn u 
Xn = e~^" = ksn ^{un + ia') sn — ia') = k 



"Ti "-^^J- ""-^ "'n 

cn ia' — cn Ur 



dn ia' + dn m„ 

^mn = -V^Sn - (3.3) 

Note that the trigonometric/hyperbolic variables, u;„ and 7„ are related to the elliptic vari- 
ables, Un and a' by the functional equation 

exp(±iiu;„ - \'^n) = \/ksn\{un± a'). (3.4) 

The choice of a' rather than a as anisotropy parameter is arbitrary because of the symmetry 
under interchange of horizontal and vertical lattice axes. Choosing a' at this point results in 
expressions for uon and 7„ which are equivalent to those of Onsager. 
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The mapping between the elhptic parameterization and the trigonometric/hyperbohc one 
was described in ref. [|11[] for the isotropic case. The formulas below are a generalization of 
this mapping. It is simpler to use the trigonometric parameterization than it is to use the 
elliptic parameterization for numerical series generation. 

The variables Un and 7„ defined above satisfy the identities given in Appendix 2 of 
Onsager's paper [Q], 

cosh7„ = {cniadnun — k diaiacnUn)/Mn (3.5) 
sinh7„ = —ik''^snia/Mn (3-6) 
— costUn = (dn ia cnun — k cn ia dnM„)/M„ (3.7) 
smuJn = k'"^ snUn/Mn, (3.8) 

with 

M„ = dniadn-u^ — A;cniacnM„. (3.9) 

Notice that it is a which appears in these formulas rather than a'. The mapping from the 
variables u„ and a to the variables u;„ and 7„ is conformal, as is seen from the formulas, also 
given by Onsager 

—dujn/dun = d'jn/da = k'"^ /Mn = 2sinh7„/ snia 
d'jn/dun = duJn/da = k''^ksniasnUn/Mn. (3.10) 

Finally, Onsager gives the functional equation 

cot\{uOn - i'^n) = (1 + fc) sc |(m„ + lo) ud |(m„ + lo). (3.11) 

To make the change of variables, we first note that (pn = corresponds to m„ = 0, and 
(f)n = corresponds to Un = 2K. Using eq. ( |3.10| ), we see that 

/7r /'4K 
duj cot a/ sinh7 . . . = y/k / du . . . , (3-12) 
■vr Jo 



which implies 

.(j) _ cot^ duji dujj_i ( Yj 1 \ ( TJ u \ l + n„=i3^»i 



•^^^ \n=l '"■/ \l<i<fc<7 / iln=l'^" 



l<i<k<j 
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The condition Ui + . . . + ooj = mod 27r, which resuhs from having performed one of the 
integrations, is assumed. In terms of Un, k and a, the quantities x„ and sinh7„ can be 
expressed as 



cot^ a 



^ — cos Un — a/(.^ — COSCUri)^ — (cot a) 



sinh7„ = cot^ a a/ — cosct;„)2 — (cot a)"^ 



with 



e = (1 + l/(fccot2 (1 + fc/cot2«)'/' = (1 + (S')-')'^' (1 + isT'f'^ 

In the isotropic case, these reduce to 



Xn = s + s ^ — COS Un — + s ^ — cosa;„)^ — 1, 
sinh7„ = ^(s + — cosa;^)^ — 1. 



(3.14) 
(3.15) 

(3.16) 



(3.17) 
(3.18) 



Here s = \fk for high-temperature and s = for low temperature, although the distinc- 

tion is irrelevant as the dependence of the integrand on s and 1/s is symmetric. Finally 



cot a 



sin \{uJi - uJk) _ 1 sinhi(7i-7fc) 
' sinh I (7i + 7fc) cot a sin | {coi + LOk) 
2(xjXfc)^/^ cot a sin 1(^74 — Uk) 

1 ^i^k 



(3.19) 
(3.20) 



When j = 1 or j = 2, the integrals can be rewritten in terms of known functions as was 
noted in the introduction for the isotropic case. When j = 3 Glasser showed that the 
integrals can be written as an integral involving the square root of a polynomial of degree 
higher than 4. Such integrals are often called hyperelliptic integrals, and are a special case 
of Abelian integrals. 

The corresponding anisotropic expressions are 



X 



(1) 



1 



X 



(2) 



(^1/2 _ ^-1/2)2 



2 CSC 2a + J (A;i/2 _ A;- 1/2)2 + 4 csc2 2a 



(3.21) 
(3.22) 



As we show in section ^, we can say more about the general form of the expansion, based 
on inspection of the long series published by Nickel 0,|lll5 as well as more recent extensions 
reported here, giving hope that there is still more regularity to be found. 
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3.2 Natural boundaries 



As mentioned in the introduction, there were observations, first in |]T^ and then in [^,TT], that 
strongly suggested the susceptibility of the Ising model is a function with a natural boundary 
unlike the free-energy or magnetization. We expand and clarify those observations here in 



the light of our new knowledge of the general anisotropic case discussed in section ^TT] and 
the additional numerical work on the isotropic limit described in section |^. 

First, we expect that as described in , the x^^^ given by the integrals (|3.1| ) multi- 
plied by the factors outside the summation in eq. ( p..9|) or ( |1.10|) will be singular at the 
symmetry points of the integrand and where the denominator factor 1 — nn=i vanishes. 
The symmetry point condition requires all uj^ to be equal and given by u;^ = ^ = Svrm'/j, 
m! = 1,2, . . .j. The vanishing of the denominator factor requires the Xn, now all equal, to 
be given by {2mm/ j), m = 1,2, . . . j. Equivalently, from the explicit formula 

(EH), 

cot2(a)(^ - cos(27rm7j)) = cos(27rm/j). (3.23) 
With cot^(a) = s'/s from ( p.9|) and ^ given by ( p.l6|) , we find ( p. 231) can be reduced to 



271771 27lTn^ 

cosh(2ir) cosh(2is:') - sinh(2ii') cos — sinh(2ir') cos — — = (3.24) 

with 771, 77i' = 1,2, ... i as discussed above. It will be noted that the left-hand side of ( p.24| ) 
is the denominator in the Onsager integral for the free-energy and thus we find the (to us) 
surprising result that the singularity of x'^^\ ^ property of the Ising model in a magnetic 
field, is intimately connected with a property in zero field. 

The full x? being a sum of x^^'\ will naively be expected to be singular at a dense set of 
points and thus have the Onsager line ( |3.24[ ) as a natural boundary. The presence of natural 
boundaries has implications for expansions about the physical singularity points s = ±1 
that are necessary to understand corrections to scaling. We briefly explore some of these 
implications but restrict ourselves for simplicity to the ferromagnetic point s = 1 in the 
isotropic model. We also make the plausible assumption that the singularity in each x^''^ 
closest to r = is the most important for determining, in expansions of the large p 
asymptotics and this considerably simplifies the discussion^. 

Let Tj = i sin 9j be that singularity in x^-'^ that is closest to the ferromagnetic r = 0; 9j 
is fixed by cos^j = (1 -|- cos<pj)/2 with 0j = 2ti/ j. Choose the branch-cut arising from this 
^Note that singularities on \s\ = 1 map to points on the imaginary axis in the complex t plane. 
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singularity to lie along the imaginary r axis and directed away from r = 0. Take t = iT to 
be a point on the branch-cut. Provided the (positive) deviation 



66 j = arcsinT — 6j (3.25) 

is not too large, the discontinuity across the cut can be estimated from the linearized singu- 
larity equations (14) in |^ and (12) in |]TT|. For j > 2 a more general result that includes 
the first order correction is 



Disc(x^^'^) = -C,2-''7sinVi[5^jSin^,/sin2 0,](j'-3)/2 



1 + 66. sm 6i < ^ H — } + 0(66f 



(3.26) 



where 



C, = 2V2{2/nY^--'y\j/2Y^'-'y'{Ul^,T{m))/T{{f - l)/2) ^ 3.7Q55r'/''2^ exp{-f /A). 

(3.27) 

The last equality is valid only in the large j limit. To obtain the discontinuity in say 
we must first sum the discontinuities in x^''\ j odd, and this we can crudely estimate by 
integrating over j, keeping only the leading exponential factors in ( 3.25|) and ( 3.26| ) and 



making a small angle approximation as well. Essentially the same formula is obtained for 
the discontinuity in x~ on summing over j even; in either case we find 

"271 /T 

'V27T/T 



1-211 /r 

Disc(^ X^^^) ~ / dj[(Tj7(v^7r) - l)/{2V~e)f'^ (3.28) 

JV2-KIT 



where the sum is over j odd or j even, and the lower limit is simply the restriction to those 
j values that contribute, while the upper limit roughly defines the limit of validity of the 
linearized approximation. The precise value of this limit is not important since the integrand 
has a maximum well below the limit. For large j only the maximum of the integrand matters 
and ( |3.28| ) reduces to 



Disc(x) ~ exp(-39. re/T^) (3.29) 

in which the numerical coefficient of 1/T^ is 7r^x^/(2(x — 1)) with x the smallest real solution 
of 21og((a; — l)/2) + l/(a; — 1) = 0. Keeping terms such as the correction term in (|3.26|) and 



the 2^ in (|3.27|) in the steepest descent analysis lead to 0(1/T) corrections to the exponent in 
( p.29| ) so we conjecture that the right hand side of ( p.29|) is exactly the leading exponential. 
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It is interesting that the discontinuity (|3.29|) is similar to that found in weak couphng field 
theory expansions, but the mechanism producing the cut here does not seem to be related 
in any way to instantons. 

The additional singularities that x has at r = 0, namely terms such as the divergent 
"scaling" |t|~'^/'^ or the "short- distance" powers of log | complication we do not 

know how to handle in any rigorous fashion^. As a consequence we will simply ignore them 
and make the simplest, yet reasonable, assumption that they make an additive contribution 
not relevant for understanding the effect of the natural boundary. Then the cut discontinuity 
( p.29| ) would imply a divergent behavior in the r expansion of x- That is to say, the coefficient 
of in the limit p — ^> oo will diverge as T{p/2) / aP^'^ with a ~ 39.76. This follows from a 
contour integral around the origin distorted to run on either side of the cut imaginary axis. 
The contribution of the cut discontinuity to the coefficient, Cp, of in the expansion is then 

Cp oc y (iT/Tf+i exp(-a/T2) ~ T{p/2)/aF/^. (3.30) 

We know little about the cut discontinuity on the circle |s| = 1 other than what we 
have deduced near s = 1 as given by ( |3.29| ). However, the fact that the amplitudes of the 
singularities of x^^'' l-^l = 1 vary dramatically with order j almost certainly implies there 
will be no cancellation of singularities in the sum of x^^^ that defines x- Furthermore, the 
variation with order means that there is no length scale at which, as one approaches the 
circle of singularities \s\ = 1, Xiso can be smooth. It is these two points taken together 
that we consider overwhelming evidence that, in the isotropic case at least, x has a natural 
boundary that is the entire |s| = 1 circle. We do not imply by this that all points on the 
circle are equally "singular". As argued above, the existence of an asymptotic expansion 
about s = 1 seems likely. A very different situation arises at a point such as s = i. While 
the singularity in x*^-'^ on the circle s = exp{i9) nearest s = 1 lies at a distance A9 = 0(l/j) 
for large j, the corresponding nearest distance from s = i is A6 = 0(l/j^). Furthermore this 
latter singularity is larger in its leading amplitude than the former by a factor roughly 



If we knew that the "short-distance" corrections Bf in (1.17) formed a convergent sequence, then a 



subtraction process similar to that described in section 6.2 could be carried out here. The numerical evidence 



from section 3.2 is that the point t = in the scaled and pole subtracted x± is no longer a branch point 



singularity and hence the final result ( 3.30 ) is justifiable and would apply directly to the scaling-amplitude 



function F± . Unfortunately, we do not have any independent information to suggest that the "short-distance" 



sum (1.15) is convergent and thus the argument leading to (3.30) is at best suggestive that one or the other 



(or both) of the "short-distance" or "scaling" sequences are asymptotic. 
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The reduction in distance and dramatic increase in amplitude suggests that an asymptotic 
expansion about s = i is not possible, but this has not been proved. 

For the anisotropic case we have not analyzed ( |3.1| ) in detail so we do not have the 
necessary amplitude information to make the same claim directly. However, we can take the 
extreme anisotropic limit of s' infinitesimal (but not 0) and find that the Onsager line ( ^.241 ) 
has now come very close to the circle |f | = 1. At this point we can connect to the work 



m 



12|. There it was observed that if x{K,K') was written as 



xiv,v') = J2Hniv)v"' (3.31) 

the Hn would be singular at a dense set of points on \v\ = 1 a.s n —>■ oo. Furthermore, if 
as in the discussion above, v' is infinitesimal (but not 0) the amplitudes of the singularities 
vary dramatically with order n and the same conclusion as in the isotropic case is reached. 
Given that the Onsager line ( |3.24| ) is a natural boundary in two extremes, it seems highly 
probable that it is also a natural boundary at all intermediate anisotropy values. 

We conclude this section by contrasting the above behavior of the susceptibility with that 
of the free-energy and magnetization which are only singular at an isolated set of points, not 
a dense set. This is precisely what one expects for a D-finite function, as in that case we 
have the following!] 

Theorem 1 Let f{x,y) = X]n>o ^"-^"(^) ^ D-finite series in y with rational coefficients. 
For n > 0, let Sn be the set of poles of Hn{x); let S = [Jn^n- Then S has only a finite 
number of accumulation points. 

This is observed in practice. The anisotropic magnetization and free-energy each have exactly 
one accumulation point , while the (non-D-finite) susceptibility appears to have an infinite 
number ||12|| . 

3.3 x^^^ is D-finite 

This remarkable result, that x^''^ (o^' equivalently x^-'^) is D-finite while x is not, follows 
from the results of Lipshitz [l^ (see also Zeilberger [Q), who gives several basic definitions 
and theorems concerning D-finiteness of series in several variables. The integrand of the 
trigonometric/hyperbolic form of x^^^ is an algebraic function of the variables cj = cos to j 
^Due to M. Bousquet-Melou, private communication. 
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and k^^'^, and is thus a D-finite function in these same variables. Then by Theorem 2.7 



of ref. |T3|, integrating over one or more variables preserves D-finiteness which implies the 
result. The term D-finite is synonymous with holonomic in much of the literature. 

Kashiwara and Kawai have shown that any Feynman diagram is holonomic, whereas 
an infinite sum of such diagrams may not be. This is just the phenomenon we observe here. 
At first glance it appears that Kashiwara's definition of holonomic differs from that used 
here, but this is not so. The point is that the definition of D-finite functions of more than 
one variable requires that the underlying system of partial differential equations be such that 
only a finite number of initial conditions are needed to specify the function. Such systems are 
called "maximally over-determined" or "holonomic" in the analysis literature. In the single 
variable case, the question of a finite number of initial conditions is clearly automatically 
satisfied. 

Motivated by the above observation, we have attempted to find linear differential equa- 
tions with polynomial coefficients in fc^/^, or equivalently, linear recurrences for the series 
coefficients, of x[^l and x[^l. With the available series of order k'^^'^^'^ and k^"^ respectively, 
calculated using the methods of numerical integration described in [3,|TT|! we have ruled out 
any such recurrences of depth 14 with coefficients of degree 15 for xl^^, and of depth 6 with 
coefficients of degree 7 for x[^l. Thus while these functions are provably D-finite, it is clear 
that the generating differential equation will be a fairly cumbersome object. 

We have also attempted to fit these series as polynomials in the complete elliptic integrals 
K and E with polynomial coefficients in k and obtained similar negative results. 

4 Correlation functions and difference equations 

In this section we give details of our more efficient method of series generation based on 
summing the correlation functions, obtained by means of nonlinear recurrences, as outlined 
in the Introduction. 

The equations we used for generation of the very long series are the ones given by 
Perk |jl5|. Here we present a slight generalization due to McCoy and Wu which keeps 



track of the separate multi-particle components. The expansions of the two-point correlation 



functions in multi-particle components are analogous to the corresponding expansions ( |1.9| ) 
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and ( |1.1CI| ) for the susceptibility 

C{M,N;X)={ ^ ' ^"=° ^ ' ' (4.1) 

y (1 - efl^ Y.n=o A2"C'(2")(M, A^) for T < T„ 

with C'(°)(M, A^) = 1 and for j > 

C^^\M,N)= (4.2) 

/ / . / TT W TT . > / TT__ 1 ---(iV^^„) 

n=l 



j! i-. 27r 




The fugacity A is associated with the number of particles, and C(M, A^) = C(M, A^; 1). 
With these definitions the quadratic partial difference equations are 

s^[C{M, N; X f - C{M, N~l; X)C{M, N + 1; A)] 

+ [C*(M, A^; Xy ~C*{M- 1, A^; X)C*{M + 1, A^; A)] = (4.3) 
s'^[C{M, N; A)2 - C{M - 1, A^; A)C(M + 1, A^; A)] 

+ [C*{M, N; Xy - C*{M, N-1; X)C*{M, A^ + 1; A)] = (4.4) 
ss'[C{M, N; X)C{M + 1, A^ + 1; A) - C{M, N + 1; X)C{M + 1, A^; A)] 

= C*{M, N; X)C*{M + 1, A^ + 1; A) - C*{M, N + 1; X)C*{M + 1, A^; A). (4.5) 

The object C*{M, N; A) is the correlation function on the dual lattice and is obtained from 
C{M,N; A) by replacing s' with s* and s with s'* . Equation ( [4.5| ) holds for all M and A^. 
When M = N = equations ([4. 3D and (|4.4|) must be replaced by 



C*(1,0;1) = VrT^-sC(0,l;l) (4.6) 



C*(0, 1; 1) = VI + - s'C{l, 0; 1). (4.7) 

We do not know of A 7^ 1 versions of these equations. 

These equations are nearly enough to determine all two-point functions completely. For 
the isotropic expansion, all that is lacking is either the high or the low temperature set 
of diagonal correlations (M = A^). From either one of these the other can be obtained 
using equation ( [4.5|) with M = N. When A 7^ 1 we have used the integral formula (|4.2| ) 
to compute the diagonal correlation functions. For the rest of this section we focus on 
the case A = 1 where two superior methods for obtaining the diagonal correlations are 
available. From the purely computational point of view the difference equations of Jimbo 
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and Miwa are almost certainly the most efficient and to be preferred. However, from 
the point of view of understanding the analytical structure of the correlations the original 



Toeplitz determinants p7|-|39| are better. We also find that for the numerical computations 
we have carried out so far the evaluation of the determinants is only a small fraction of the 
total project time so efficiency is not yet an issue. 

We will restrict ourselves in the following to the isotropic lattice. In that case and for 
> 0, the diagonal correlation C{N,N) is the determinant of an x Toeplitz matrix 



with elements a 



that are the integrals |39(| 



an = i2n)-^ / d9{s - exp {-2ie)/s) exp (-2i^n)/ VT^Tsh^. (4.8) 
Jo 

These integrals apply both above and below and furthermore since r ^ — r corresponds 
to s —>■ 1/s, one can establish from the integral in ( [4.8| ) the relations a-n-i(T) = — a.n(— t)- 
Explicit formulae for a„ for small n can be given in terms of elliptic integrals E and K but 
these expressions are not particularly enlightening and are not needed here. Rather we need 
the series expansions of (|4.8| ), either in s, 1/s or r depending on the application. 

The series expansions in s and 1/s both for a„ and C{N, N) are completely straightfor- 
ward with computer packages such as Maple that automatically handle the multiple precision 
arithmetic required. Furthermore, these same packages can be set to treat the C{M,N) in 
the recursion formulae (|4.3| )- (|4.5| ) as series and thus very little programming is necessary to 
generate the susceptibility. Admittedly, some steps need to be taken to conserve time and/or 
memory resources but this is very hardware dependent and will not be described here. What 
is worth noting, however, is that for a series to order the recursion formulae require 0(A^^) 
multiplications of series of length A^ and thus in a naive implementation, 0(A^^) multipli- 
cations. Since the word length grows linearly with A^ the algorithm has complexity of at 
most 0(A^^). There are more efficient ways to multiply long series and numbers with a large 
number of digits |^ but we have not found it necessary to explore these options. 

Timing proportional to A^^ is observed in practice, in our implementation of the recur- 
sion in Maple. We have generated high-temperature series of order 323 and low-temperature 
series of order 646 for the isotropic susceptibility. The entire calculation took 123 hours 
on a 500MHz DEC Alpha with 21164 processor running Maple V version 5.1. We have 
also obtained shorter anisotropic series in this way (either the nearest off-diagonal correla- 
tion functions, or additional assumptions are needed). As lattice anisotropy is a marginal 
operator, we hope that an extension of this calculation will be very revealing. 
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The series in r is most easily obtained by expressing the Toephtz element integral ( |4.8| ) 
in terms of hypergeometric functions. To show this connection we start by writing 



Tca- 



,/2 = Vl + tM, - tA, 



(4.9) 



in terms of the real integrals 



A.. 



Ar 



tt/2 



7r/2 



de sm{e) sin(z/^)/\/r2 + sin2^, 
de cos{e) cos{ue)/\/T^ + sin^e, 



(4.10) 



where u = 2n + 1. The required symmetry a{—v,T) = —a{v,—T) is explicit in eqs. 



and (^4.101) and by standard integration by parts manipulation one can show that the cosine 
integral satisfies the differential equation 



(1 + T^){d/dT)T{d/dT)A^ - u^tA^ = 
while the sine integral can be expressed as the derivative 

As = -{T/u){d/dT)A,. 



(4.11) 



(4.12) 



Furthermore, a direct evaluation of the integral in ( [4.10 ) for small r yields A^ = —1^ + o(r) 
with 



4 = log(|r|/4) + ^(z//2)/2 + V^(-z//2)/2 - V^(l/2), 



(4.13) 



and this initial condition together with ( |4.11| ) completely determines Ac- Since (|4.11|) can 
be recognized as the hypergeometric differential equation in the variable z = — r^, we can 
immediately write the solution as 



A = - i.F{u/2, -vjl- 1; -r^) + ^(zy/2),(-z//2),/(A:!)^(-r2)'= 

X (^/.(A; + 1) - + i){v 12)12 - ^(k + V 12)12 + ^(-z//2)/2 - i)(k - v/2)/2). 

(4.14) 
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On Taylor expansion we now obtain 



7ra„/2 = + 4r + [z/(4 - 1/2) + VV2 + y\L - l)rV4 (4.15) 
+ [z/3(4 - 5/4) - 1/ - 2z/- VVl6 + K(4 - 3/2) - 4z/2(4 - l/2)]rV64 + 0(r6) 

which can be extended as required. 

From the expression ( 4.15 ) one can conclude that C{N,N), as an N x N determinant, 



will contain logarithmic terms r^(log|r|)P with q > p and p < N — barring cancellation. 
However, there is cancellation and the key qualitative features of the cancellation can be 
deduced simply by using an alternative representation for the determinant. In particular, by 
systematically subtracting rows and columns one can show that an equivalent determinant 
has matrix elements aij which are of the same integral form as the a„ = aij in eq. ( [4.8|) except 
for the replacement of the exponential exp{—2i6n) by the product exp(— z0n)(2 sin(^^))*+-'"^. 
The new matrix is no longer of Toeplitz form but for our purposes here it is the better 
representation because of the powers of sin(6') which have the effect of shifting the log|r| 
singularity of the integrals to higher order in r. Indeed one can show from the new integrals 
that the leading singular behavior of each matrix element aij is proportional to r*^-'"^ log |r| 
and this is sufficient to show that in the logarithmic terms r''(log |r|)^ in C{N, N) one must 
have q > p"^. We have not attempted to pursue this argument to deduce C{N, N) analytically 
but rather have resorted to a numerical small evaluation of C {N, N) using ( [4. 151) and then 
fitting to obtain formulae valid for general N. Our results for the leading logarithm term are 
summarized by the expression 



^/^C{N,N,T) = C{N,N,T = 0)J2^''i^og\T\+CNY{Nr/Ay' (4.16) 

p=0 

p-i 

{Y[{N-^ - k-^f + (1 + 2{N^ - p^)y/8 + 0(r=^)} 

k=l 

in which we have used C to denote the discrete logarithm, i.e. 

£jv = ^{N + l)/2 + ^(iV)/2 - ^(1) - log(4). (4.17) 

The product factor in the first braces of ( |4.16|) is to be understood as unity for p < 2 and in 
the final brace pair the coefficients of are polynomials in N of degree < q. Furthermore, 
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with the ^/s factor extracted exphcitly as in ( [4.16|) these coefficients of r'^ vanish if g = 2/c + 1 
with k < p. The critical correlation factor in ( |4.16|) is 

N 

C{N,N,T = 0) = l[T\n)/{T{n + l/2)T{n-l/2)), (4.18) 

n=l 

and approaches A/N^^^ as N ^ oo where 0] log(y4) = 3C'(— 1) + log(2)/12. We have 
used ( [4.16|) , extended or truncated to some order in r, as input to the quadratic recursion 
formulae to generate all correlation products y/sC (m, n) as series in r within an octant 
■m>n>0, m + n< 2N + 1. While most of our results are numerical, they are consistent 
with the assumption that the structure observed on the diagonal remains true on the octant. 
That is, if we define n = fiN, < /i < 1 and set m + n = 2N (even shell) or m + n = 2N + 1 
(odd shell) then on these even/odd shells for large A^, the correlations are of the form 

oo 

V~sC,/o{f^) = iV-V4 ^(log |r| + C^YiNr/AfA'^^^ifi) (4.19) 

where the A^^^^{jj) are still Taylor series in r as in ( [4.16| ) but with coefficients that are now 
(possibly asymptotic) Laurent series in 1/A^ rather than polynomial in N. The highest power 
of multiplying t'^ remains A^*^. 

Considerable care must be exercised in the numerical work since the recursion formulae 
are unstable. A toy recursion relation of structure similar to the ones we use in the Ising 
study illustrates this nicely. Let 



Cm+l — (2(Cm) ~ ) /^m-l (4.20) 

which is to be applied to all possible n and iterated forward in m. The recursion ( [4.20|) has 
as a solution a constant, say M, but is susceptible to a steady state growth of errors so that 



^M + e(-l)V 



(4.21) 



is also a possible solution. On substituting ( 4.21 ) into ( [4.20 ) one finds the growth constant 
a must satisfy — 6a + 1 = so that 



a = 3 + VS, logio(«) = 0-765. 



(4.22) 
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This is very close to what we find in our numerical work and implies that if we want some 
number of digits D that are accurate at the outer edge of the octant on a shell specified 
by + m = 2N + 1 we need to start with digits Dq D + 1.53A^. In practice we have 
worked as high as = 146 with Dq = 380 using the automatic multiple precision facility 
of Maple. In section we use these results to calculate the "short-distance" (including 
analytic background) terms in the susceptibility. 

5 Conjectured short-distance structure 

In this section we state some conjectures for the short distance behavior of the A 7^ 1 model 
introduced in section ^. We arrived at these conjectures by inspection of series obtained from 
a combination of the integrals (|4.2| ) and the difference equations (f4.3|) - ([4.5|) . Our interest 
in this model is motivated by several considerations. Firstly, it enables us to see how the 
analytic structure of x evolves as successive contributions x*-^"''' or x^'^^'^^^ are added. Secondly, 
the correlation function C{M, N; A) for small values of M and may be required as initial 
conditions for certain series generation algorithms. Thirdly, we hope that the presence of 
an additional parameter which can be varied will provide some insight into the Ising model 
itself. Finally, the deformations of the elliptic functions that appear in our conjectures may 
be of intrinsic mathematical interest. 

Following the lead of ref. we make a change of variable from the modulus, k, to the 
nome q. Examination of the j'-particle contributions to the isotropic susceptibility and two- 
point functions reveals that there is much regular structure. We arrive at exact conjectures 
for (7(^^0,0), C^^\l,0), C^^\l,l), C'(^)(2,0) and C^^\2,l) as functions of j and q. These 
provide the first terms in the short distance expansion of the susceptibility. 

5.1 g-series in the Ising susceptibility 

Although more complete results have been obtained for the correlation functions, we will 
demonstrate the method by which we derived our conjectures using the susceptibility as an 
example. To make our observations more concrete, we reproduce tables of series coefficients 
for (1 — k'^y^^xlso- an example of how to interpret the table, we read (1 — k'^Y^^X^gl = 
8(A;9/V29 + 4A;V2^2 + i6A;W2/2i3 + 4^k'' /2^^ + 20k''>/y2^' + Mk''/2'^ + ...). 

Inspecting Table |I|, we make the conjecture that the series can be sensibly decomposed 
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Table 1: Coefficients of the series (1 — k^y^'*'x[il/2^ . The numbers in the tables are the 
coefficients of {Vk/2)'^ starting at n = in the upper left corner and reading down and to 
the right. 
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1 


126 


8760 


AAA ^7 A /~\ 

444740 


-I A £^ A £^ 

18429842 


661181352 


j 


— 8 






























4 


584 


46440 


2666700 




























1 


160 


13839 


858704 


42821009 


1823591632 


3 


= 9 








4 


708 


67252 


4553260 


















44 


8552 












4 


728 


70976 


4924124 






1 


198 


20888 


1560492 


92610504 


4644898080 


3 


= 10 






























4 


868 


99812 


8087916 




























1 


240 


30359 


2692864 


188045229 


11006289872 


3 


= 11 











4 


1064 


-^ A C\ A C\ A 

148424 
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4 


1044 


1 A r\r\r\ 

143020 


13686020 






1 


286 


42756 


4447860 


361695338 


24490780096 


3 


= 12 

































4 


1256 


205352 




























1 


336 


58631 


7076256 


663817077 


51575531568 


3 


= 13 


u 








4 


1/1/1/1 

1444 


270020 





































4 


1464 


277424 






1 


390 


78584 


10898664 


1169440708 


103475590040 


3 


= 14 

































4 


1668 


358612 

























Table 2: Table continued. 
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as 



2^{^/2f[{l + co,iA;2 + co,2A;" + . . . ) 

+ (v^/2)^'(4 + k + ci^2k^ + Ci^^P + ...) 

+ {Vk/2f^{c2,o + C2,ik + C2,2k^ + ...) + •••] (5.1) 

In fact, we observe that as j tends to larger and larger values, the first row of coefficients 
tends towards the expansion in k of 2^q^^^^/63{0,q)!^ The first correction comes in at order 
Hence we make the change of variable from k to q in x[^l and divide the result by 
2^q^^/'^/9s{0, q). For all j we obtain a series of the form 1 + Aq^ ^ + cq^^^^^ ^ + . . . The sequence 
of terms starting at order q^^^ appears again to be fitted by a recognizable function of g, at 
least until contributions appear at order q^^^^ or g^-'/^. Subtracting this assumed product 
form yields a series whose first correction enters at order q^^/^. The coefficients of the terms 
of orders lying between g^-'/^ and q^^/'^ are not independent of j as was the case previously, 
but vary linearly with j. The constant part has a product form, and the j-dependent part 
may as well be we do not have sufficiently many terms to make a firm conjecture. Likewise 
the terms between q^^^^ and q'^^f^ depend quadratically on j, and the correction at q^^^^ 
appears to vary as the fourth power of j. 

From the long series we have produced, we have been able to conjecture that 



CO 



1 + g 



ri-l/2 



n=l 




^3(0, g) 

oo 

11(1 - q'^-ni + q'^-'^r + ij + iW^VOliO, q'/'^ 

n=l 

ff\q) + (j + 1)(J + 2)/f (g)) + O [q'^"] 



(5.2) 



where 

ff^ = 1 + 6gi/2 + 26g + 17g^/^ - 81g^ - 55g^/^ + 285g^ . . . (5.3) 

/f ) = 1 + gl/2 _ 5g _ 4^3/2 ^ ^5^2 ^ ^Q^5/2 _ 3g^3 ^ _ _ _ ^5 4^ 

These are consistent with the expressions for the correlation functions in the following section. 
That is to say, summing the correlation functions given in the following subsection gives terms 



that agree, as far as they should, with the above expression. Similarly, summing (|5.2|) over j 



gives a series that agrees to the appropriate (low) order with the known expansion for Xiso. 

^This may be checked by reverting the product form given in (|]^) for k = k{q) and using the identity 
^3(0; 9) = ^/SK/tt to obtain the expansion of ^3(0, q) in k. 
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5.2 g-series in the two-point functions 



In this subsection we write down some conjectured results for the short-distance correlation 
functions which we have derived empirically. For these cases, unlike the susceptibility, the 
general correction term is apparent from the series and we are able to formulate exact 
conjectures. Let us define operators, $o a-nd which convert power series in z to power 
series in q according to the rules 



n=0 n=0 



(5.5) 
(5.6) 



n=0 



n=0 



Then we conjecture that 
2-^(1 -A;')^^^^) (0,0) 

2-^(1 -A;')^/^^^^) (1,0) 
2-^(1 -A;')^/^C'(^)(1,1) 

2-^(1 -A;')^^^^) (2,0) 
16 



1 z^l-z'') 



^3(0, g) '(l + ^2)i+i' 

1 / 2g^/^g3(0,g'/') \ '''' ^ z\\-z) 

^3(0, g)V ^2(0,gl/2) ) l(l + ^2),+i: 

Z^'+^(l - Z^) 



2(J + 1) 



^2(0,g)^i(0,g) (1 + ^2)^+2 ' 
1 / 2gVsg3(0,gV2) y/2 



(5.7) 

(5.8) 
(5.9) 



gV4^3(0, g) V ^2(0, gV2) ; ^ (1 + ^2)j+l(l _ ^2) 



,1/4 _ 



^3(0,g)^|(0,gV2)^ rfgV4 
^2(0, g) 



$0 



^J+2 



(1 + ^2)^+2(1 _ ^2) 

8 



^3(0, g) 



1/4 _ 



2^32(0, ?)^2(0, g^) ^3(0, ?)^|(0, ?V2)^2(0, ^4) - ^^1/4 



X $r 



^j+2 



(1 + Z2)^-+2(1 _ ^2) ' 



and 



2-^(1 -A;2)V4(70)(2,i) = 



^3(0, g) 



8g 



^1(0, gV2) ^^1/4 



$1 



(1 + ^2)2(1 _^2)_ 

(5.10) 



Z^+^(1 + Z) 



(l + ^2),-+l(l _^2) 



+ ^1 



^^-+^(1+^) 



i-z^ 



$1 



X 



0,5/16 ^2(0,?^/^) 



(l + ^2y+l(l _^2) 

^(1 + ^) 



$1- 



(5.11) 



1-z^ 

As noted above, these are necessary, but not sufficient, for the generation of the j- 
particle contributions to the correlation functions, being some of the initial conditions for 
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the recurrences. 



6 Scaling form of the susceptibility. 

The main result in this section is a conjecture completely specifying the analytic structure 
of the susceptibility x of the isotropic Ising model in the vicinity of the critical point both in 
the neighborhood of the ferromagnetic point s = 1 and the anti- ferromagnetic point s = — 1. 
The conjecture, contained in eqs. ( p. . 1 5| - [r!T^ ) is based on what we believe is overwhelming 
numerical evidence that is obtained by disentangling the "short- distance" and "scaling" 
parts of X in a manner described below. 



In section |6.1| we give the assumptions and numerical procedures we use to derive the 



"short-distance" part of x ( |1.15 ) with the coefficients listed in the Appendix. Then, in 



section |6.2| we describe the "short-distance" subtraction and analysis on which the behavior 
of the scaling-amplitude function F± shown in ( |1.18 ) is based. Finally, in section |6.3| we 



outline our fitting programs to determine the coefficients in the functions F±. Since there 
are no confluent singularities in F± whose amphtudes need to be found, the fltting procedure 
is very well-conditioned and the coefficients in the Appendix are as determined to an accuracy 
of up to 20 digits. 

6.1 "Short-distance" term. 

That the "short-distance" contribution to x can be obtained from numerical values of C(m, n) 
for small |m| and \n\ relies on certain assumptions about the behavior of the expansion 
coefficients in eq. ( [4.19|) . In particular, we assume that ( ^4.19|) remains valid up to of 



the order 1/r where it can, in principle, be matched term by term to a large distance 
expansion that properly describes the roughly exponential exp(— A^r) decay of correlations 
as A^ ^ oo. Explicit matching formed the basis of the previous calculations of terms in the 
"short-distance" x (cf- 01) ^^^^^ becomes extremely cumbersome at higher order. Our 
ability to go to high order here rests on the fact that we dispense with such explicit matching 
and rely instead on power counting to uniquely identify those terms that contribute to the 
"scaling" and the "short-distance" parts of x separately. We believe this is tantamount to 
the scaling argument that in the critical region there is a single length scale proportional to 
l/r'^ with V = 1 and thus, in a way to be made more precise below, we can deduce the power 
law of the large distance contribution of any set of terms varying as A^^ by simply replacing 
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A^^ with Terms whose variation is as a fractional power of r (with possibly logarithmic 
multipliers) are discarded as assumed contributions to the "scaling" part of x- Terms whose 
variation is predicted to be an integer power of r (with possibly logarithmic multipliers) are 
assumed to be part of the "short-distance" x ^^id are treated more carefully. 

To make the argument and assumptions more explicit we begin with some definitions that 
will be useful also for the subsequent analysis. Let the two dimensional sum (|1.8|) defining 
X be reduced to a one dimensional sum by combining the contributions from all sites on 
the even and odd squares \m\ + \n\ = 2N, 2N + 1 and then further combining these into 
sum and difference combinations which are necessary for separating the ferromagnetic and 
anti-ferromagnetic contributions. That is we write 

V~SC^, = Y^d^y = V^5^(Ce(/x) ± a(/^)) (6.1) 
k 

where the first sum simply defines the coefficients while the second is the actual lattice 
sum over the octant correlations defined by ([4.19|) and extended by symmetry over the full 
square. The coefficients A^j^'f^ in the expansion 



CP, = i2^Nl\^og\r\+Cr,y (6.2) 



have the large N asymptotic expansions 

p=0 p=0 

as assumed in ( [4. 191 ) based on numerical evidence. The ferromagnetic and anti-ferromagnetic 
cases have been treated separately in ( |6.3|) to emphasize that shell subtraction in (|6.1|) reduces 
the leading power of by unity. The upper limit L on the logarithmic powers in (|6.2| ) 
depends on k as discussed in connection with eq. ([4.16|) but its precise value is not needed 
in the following discussion. We also introduce the partial sums 

N L L 

n=0 1=0 1=0 

where the R^^j'f^ in the second equality in (|6.4|) are numerical coefficients that are directly 



generated by the quadratic recursion relations and the subsequent lattice summations. In 
the final equality in (|6.4| ) these expansion coefficients have been formally split into two terms. 
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This formal separation is to be understood in an asymptotic sense with (3^'''^ a constant that is 
the N independent part of -Rj^'^"* as — * cxd and which we can understand as an "integration 
constant". The remainder 6Rfl'^\ we will show in the case that / = L, has an asymptotic 
expansion like the A^'^'* in ( |6.2| ) except for an extra factor of from the summation. For 
can be expressed as sums of such expansions with multiplying logarithmic 

factors Cn- 

We now come to our key assumption. The partial sums in ( |6.4|) in most cases diverge as 
N oo because of the presence of large fractional powers of in the 6R^lj'^\ However, we 
assume that if the formal matching of the short and large distance expansions for C (m, n) had 
been carried out, these partial sums would in fact converge and furthermore the contribution 
of each power of term could be estimated by the replacement A^ — > 1/r. Since all the powers 
of N in (5/25i'f 

are fractional, the result is that these terms contribute only to the "scaling" 
part of X- Thus the "short-distance" part of x comes entirely from the "integration constant" 
term jS^'''^ in (|6.4|) , and on explicitly reinserting the factors of ^/s and we get 

L oo 

^f/af = EE^±'''^'(l°gl^l)Vv^= EE^l/a^^(l°g|^l)' (6.5) 

k 1=0 p=0 q 



at the ferromagnetic and anti-ferromagnetic points. The last equality in ( |6.5| ) is eq. ( p..l5[ ) 



and follows simply by expanding ^/s in a series in r. The rest of this section comprises a 
discussion of the numerical scheme we use to isolate the "integration constants" efficiently. 

A technical problem arises in that one must somehow isolate the constant 13^^^^ from 
an A^-dependent sequence R^}^^ and this can be a very unstable procedure if the R^}^^ are 
combinations involving logarithms. Fortunately there are no logarithms in R^^^^ and one 
can iteratively remove the logarithmic factors in the remaining R^^^ by working in sequence 
from / = L, L — 1, ... to / = 0. To understand precisely how to do the subtractions we first 

il k) 

relate the unknown A^— dependent structure of the Rjs^J to the known, or rather assumed, 
simple structure of the ^j^^''. From the definition of S^^^ as a partial sum it follows that 
5*^1 — 5*^11^ = C"^! and by comparing coefficients in eqs. (|6.4| ) and ( |6.2| ) we get 

<^ - = E (7) (6.6) 

m=l ^ ^ 

In the case that / = L there are no logarithmic factors on the right hand side of (|6.6| ) and 
the equation is just a discrete first order differential equation in A^ whose solution is an 
integration constant additive to an asymptotic series of the same form as in (|6.3|) except for 
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an extra power of A^. The easiest way to determine the "integration constant" (3^ is to 
fit it together with unknown coefficients defining the asymptotic series part to a sequence of 
R^j^f. This is a stable numerical procedure, although if k is large, (3^''^^ is sub-dominant 
to many much larger terms and it is essential to calculate in high precision. For example, 
we might start the recursive calculation of R^^f^ as described in the previous section with 
Do = 380 digits, end with D = 155 digits at = 146, and continue with this number 
of digits in a 71 X 71 matrix inversion to obtain the (3^'''^ of interest while discarding the 
remaining 70 coefficients of the asymptotic series! The results are slightly more accurate 
than the final answers that are given in the Appendix. 

For / 7^ L the presence of logarithmic factors on the right-hand side of (|6.6|) makes fitting 
to R^lj'f^ impractical. Instead we define subtracted functions 



Fiif = + E ( 7 ^^^^'^ - p±''^)(-^Nr~' (6-7) 

which we use as replacements for -R^^'^'' in all fitting procedures. That the fitting functions 
^N^^ will give the same (3^''^^ is obvious since the added terms in ( |6.7| ) all contain logarithmic 
factors and are thus not A^- independent. That the F^^'''^ are logarithm free follows from 



m=l+l ^ ^ 

and an argument by induction starting from / = L. Eq. ( |6.^ ) can be verified using the 
definition ( |6.7|) and the subtraction eq. ( |6.6|) . 

To summarize, we obtain (3^^^^ numerically by using ( |6.7| ) iteratively, starting with / = 
L — 1, to generate F^''^ using procedures similar to that described above for the initial F^^'^ 
which are just the R^^f^. The results are given in the Appendix. 

Our arguments above are at best plausible and to provide a rigorous proof justifying our 
numerical procedure we believe one would have to do three things. First, one would have to 
show the C (m, n) expansion has the form ( |4.19| ) we deduced on numerical grounds. Second, 



one would have to show there exists a corresponding asymptotic expansion valid at large A^. 
Third, one must show both expansions have a sufficiently large domain of validity that the 
matching we assumed could in fact be carried out to arbitrarily high order. On the other 
hand, we have numerical evidence for the validity of our procedure as described in the next 
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section. The "short-distance" terms we have calculated, when used in a series subtraction 
process, leave a residual in high order series coefficients that is consistent with the complete 
elimination of all "short-distance" terms, both singular and analytic, to the 0(r^'^) we have 
worked. While the main intent of section is actually quite different, it is highly unlikely 
that the cancellations necessary to yield the expansion of the scaling- amplitude function 
F± of eq. ( [1.18| ) in pure integer powers of r would have occurred had there been an error, 
numerical or otherwise, in the results of this section. 

We conclude this section with a toy model example to illustrate the possible convergence 
properties of the short- distance expansion. The main impediment to extending the calcula- 
tion described above is that we do not have simple analytical expressions for the r expansion 
of the C{M,N) in general. An exception is on the diagonal where we know, cf. eq. ( |4.16| ), 
that the coefficient of r^^ilog |r|)P is 



^(P) ^ 4P-p'c{N, N,T = 0)NP - N^k 



2\p-k 



(6.9) 



k=l 



and a simple expression for the asymptotic expansion of C{N, N,t = 0) at large N can be 
found in Au-Yang and Perk [^. We can define the diagonal partial sums S^^^ = Yln=Q^^^ 
and, as in the analysis described above, ask for the short- distance coefficient lif^ which is 



the N independent part of S'^'* in the limit ^ oo. This term 6^^^^ could be viewed as a 
partial contribution to the short- distance terms of interest but since we do not know what 
cancellations will occur when we include all C(M, A^), we prefer to consider it only as a toy 
result that is suggestive for the convergence of the short- distance terms with order p. 



dp) 



The term 6^^^g is also equal to the e independent part of Yl'^=Q^n' e~'''^ in the limit 
e — * 0. This latter expression is more convenient since we can add to it any term such as 
1/(1 — e~^)P^+^/^~'^ for integer k without contributing to any e independent term. If we 
expand such terms in series in e"*^, we obtain as an equivalent sum 

K 



{p). 



E 

n=0 



k=0 



(6.10) 



Now choose the Qk in (6.10) such that the divergent terms in the asymptotic n expansion 
of Cn^ cancel the divergent terms in the Gamma functions. Then the n sum becomes 
convergent even with e = 0. With the cancellation extended to include also some slowly 
decaying terms in n, thus requiring K > p'^, one obtains the explicit formula 



diag 



n=0 



(jip) 



K 

E 

fc=0 



gkT{n + p^ + ?,/A~k)/n\ 



(6.11) 
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which is very convenient for numerical work. We find the fc^j^^g calculated from ( |6.11|) for 
p = 1, 2, . . . , 20 are 

-1.30 X 10-\ 4.29 X 10-^ 1.58 x 10'^, -9.12 x 10"^ 

-6.31 X 10~", 8.31 X 10"^^ 2.38 x 10"^^ -3.73 x 10"^^ 

-2.08 X 10"^^ 1.18 X 10"^^ 3.61 x 10"^^ -1.78 x 10~^^ (6.12) 

-6.67 X 10"", 5.80 X 10"^ 5.15 x 10"\ -1.44 x 10+^ 

-5.34 X 10+^^ 8.13 X 10+2^ 2.06 x 10+^^ -2.70 x 10+^^. 



The sequence in ( |6.12 ) clearly shows the asymptotic nature of the toy expansion. Fur- 
thermore the magnitude of the terms is in semi-quantitative agreement with r(p^/2)/a^^/^ 



from eq. ( |3.30 ). Thus while we cannot conclude that this will be how the short-distance 



susceptibility terms -Bf/af in ( p..l5| ) will behave, it is at least encouraging to note that there 
is no evidence for any behavior more singular than that predicted by our natural boundary 
analysis. 

6.2 Series "proof" of F± behavior 

Any numerical analysis program to deduce functions such as the scaling-amplitude functions 
F± from their series expansions is in essence a fitting routine and always presupposes a 
knowledge of the analytic structure of the result. Although we have very little exact knowl- 
edge of F± we do know that the multiplier of |t|~^/^ in x^"^^ (cf- eq. ( |1.13| )) contains log|r| 



terms and there is numerical evidence that this is true of the x*^"-* for n > 2 also. On 



the other hand, Gartenhaus and McCullough p2[ found that the series for x were consistent 
with the absence of logarithmic corrections in through order and this was subsequently 



confirmed for F_ as well |Tl[]. Scaling arguments on the question of logarithmic corrections 
are necessarily inconclusive because of the lack of information on amplitudes, some of which 
may vanish. We know of no scaling argument that either definitely requires the presence of 
logarithmic terms or can definitely exclude them. 

The natural boundary arguments given in section preclude the possibility that x has 
a convergent rather than asymptotic expansion about r = 0, but we know of no analytical 
argument that shows whether this applies to both the "scaling" and "short-distance" terms 
or just to one or the other .0 Neither could we distinguish the convergent from the asymptotic 
simple ratio analysis of the available terms in the expansion of ^/sF± in the variable (see appendix) 
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expansions from the data we have. Ahhough there are singularities in x dense on the hne 
—i<T<i, the singularities close to r = are extremely weak and we have not detected 
any singularity closer to the ferromagnetic r = point than that arising out of x^^^ at 
9(r) = \/7/4 ^ 0.66. Thus for all practical (numerical) purposes we can ignore the possible 
asymptotic nature of the "scaling" terms. 

The possibility that the F± might be expanded in a series in integer powers of r can 
be confirmed numerically without evaluating any of the coefficients in the expansion. The 
trick is to generate the series for the susceptibilities divided by the factor (1 — k'^Y^^. Any 
term in or F_ that is not a pure integer power law t^, p > 1, will necessarily contribute 
to the high order coefficients of the series. Of course, the leading two terms 1 + t/2 in F-t, 
which occur as poles proportional to and 1/r in the scaled x±j will also contribute, 
but the amplitude of these terms is known and so their contribution can be subtracted. 
Similarly, given the accurate "short-distance" amplitudes in the Appendix, the contribution 
of all "short-distance" terms can be similarly eliminated through order r^^. We find that the 
remaining high order scries coefficients are plausibly consistent with the "short-distance" 
0{r^^) that has not been subtracted and thus that there is no numerical evidence for any 
powers of r other than pure integers in the scaling- amplitude functions F±. The rest of this 
section gives details of the analysis that is the basis for this conclusion while the following 
section reports on our procedures for estimating the coefficients in the r expansions of F±. 

As outlined above, our search for possible terms other than those with pure integer powers 
of r in the F± involves the observation of the high order series terms in the scaled and pole 
subtracted susceptibility functions 

Ax+ = /3-'x+/(l - - (2K,V2f/'C^2-'/y(l - s)^ (6.13) 



m=l 

applicable for T > Tc and T <Tc respectively. Since the procedures for T > Tc and T < Tc 
are not different in principle, we will restrict the discussion below to the T < Tc case only. 

To implement the "short-distance" subtraction we can restrict ourselves to determin- 
ing bn, the contribution to the coefficient of arising from the — 1 singularity in 

is indicative of a convergent series with radius of convergence « 1.2. This impHes a conjugate pair of 
singularities at r « l.li. If however the series is asymptotic, this observation will fail to hold. 
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(1 - s-^)-^/^(log(-r))P/v^.0 The contribution from r'?(l - s-^)-^'^{\og{-T)f / follows 
trivially from 6„ by q repetitions of the derivative operation Drbn = (&n,+i — 6„-i)/2 because 
of the simple form r = — s)/2. Note also that the contribution from the 1/s = — 1 
singularity is identical except for an overall sign (—1)" and thus simply requires that in the 
end we set all odd power amplitudes to zero and double the even ones. To determine 6„ is a 
standard exercise in complex variable contour integration; we deform the contour in an 
integral with integrand s"'"'"^/(27ri)(l — s~^)~^/'^{\og{—T)Y/^/s to surround the branch-cut 
1 < s^^ < oo and in a final step set s^^ = expo;. The result is 

POO 

bn = 7^^^ / (ia;exp(— nx)(2sinh2a;)~^''^3{exp(i7r/4)(logsinha; — ztt)^} (6.14) 
Jo 

and although the integral (|6.14| ) cannot in general be done in closed form, expansions valid 
asymptotically for large n are easy to generate with computer algebra packages such as 
Maple. Because our series are particularly long these asymptotic expansions are essentially 
exact and are also very convenient for the subsequent differentiations. 

It is also necessary to subtract and/or smooth out the contributions from complex sin- 
gularities on the circle |s~^| = 1. The only significant ones for the present calculation are 
those at = ±i and ± exp (±z7r/3). The subtractions of the leading contributions from 
the latter singularities arising from x^^'' can be obtained directly from eq. (28) in [|lT|. The 
result is 



+ 1536/(50057r) (4 V2/(2m) !) [3^/Sin(2m7r/3 + 7r/4){(-13/2)2„j 
+ 13/4(-15/2)2^ + 1265/1632(-17/2)2„ - 5365/384(-19/2)2^} 
- 3-1/^ cos(2m7r/3 + 7r/4){5/2(-15/2)2„ + 75/8(-17/2)2,„ 
+53321/20672(-19/2)2™}] (6.15) 

where (z)„ = r{z + n)/r{z) is a Pochhammer symbol. 

We find that the further smoothing to give the remainder 



i?; = n^^Dtn-'/^n^Dtrn^K- (6.16) 

^"^The denominator ^/s factor here and in the foUowing is included for the convenience of allowing us to 
use the "short-distance" amplitudes exactly as tabulated in the Appendix. It also simplifies the final formula 

(Hi. 
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is more than adequate. Here as in [|rT| Dagn = {gn-i + 5'n+i)/2 suppresses the contributions 
from = ±i while D^gn = {gn-2 + gn + 5'n+2)/3 reduces what has not been subtracted by 
( |6.15| ). Note that the smoothing in eq. ( |6.16| ) scales up the amplitude contributions from the 
singularities at = ±1 by a factor n^^. 

We find that as we include higher and higher orders of the "short-distance" terms in 
the subtraction process the remainder amplitudes ( |6.16|) decrease in a smooth fashion. The 
amplitude R~ we obtain after having subtracted all terms through 0(r^^) listed in the 
Appendix is, a.t n = 600, ~ —4.1 x 10^ compared to ~ 1.5 x 10'^'^ one gets without pole 
subtraction in ( p.l3| ). The residual is very reasonably the amplitude we would expect from 
the 0(r^^) terms and this has been confirmed by extending the "short-distance" amplitude 
sequence in the Appendix by two terms using a crude Fade analysis. The result of this 
additional subtraction is to reduce R~ at n = 600 to ~ (—1 to 1) x 10®. Analysis of the 
high temperature series leads to a similar conclusion and also that the anti-ferromagnetic 
"short-distance" terms in the Appendix represent the susceptibility at this point completely 
with nothing left out. 

Although our analysis does not "prove" the absence of powers of r other than pure in- 
tegers, we can put very stringent bounds on the amplitudes of any possible singular terms. 
To make this comparison concrete, we suppose either of F± contains the singular term 
ApT'^ log |r I and for simplicity assume p integer. The contribution of this term to the am- 
plitude of the coefficient of s" for T > or 1/s" for T < in the high/low temperature 
series expansion of (|6.13| ) will be about ApTlp — l)/n^ relative to the pole contribution and 
this is to be compared to the observed amplitude we obtain after subtracting or smoothing 
away the known singularities as best we can. Our current results are the amplitude bounds 



\Ap\ < 10-^^3007r(p- 1), T > T„ (6.17) 
\Ap\ < 10-^^6007r(p - 1), T < Tc, (6.18) 

and these bounds essentially exclude any singularity with reasonable amplitude, scaling as 
T^, for all p less than about 15. 

6.3 F± coefficient analysis 

The task of determining F± numerically is enormously simplified by the a priori knowledge — 
strictly speaking a conjecture based on the numerical work of the last section — that F± has 



37 



an expansion in integral powers of r near r = 0. The absence of any confluent terms means 
that many different analyses can be used efficiently and we report here on two independent 
calculations that give essentially identical results, thus again confirming the above conjecture. 
We have not seriously attempted to optimize our analysis to give the most accurate numerical 
values possible and thus if it ever becomes necessary one could almost certainly improve on 
our coefficients as given in the Appendix. 

One particular feature of the F± expansion is worth noting here. When the functions are 
scaled by ^/s the resulting series are numerically consistent with series in even powers of r 
only. We have also noted a similar simplifying role played by ^/s in the "short-distance" terms 
which ultimately trace back to the Toeplitz determinant giving the diagonal correlations 
C{N,N) and to the quadratic recursion relations for general C{M,N).f^ This result is not 
entirely unexpected. The singular part of the free energy in zero field is an even function 
of r and the magnetization = (1 — s~^y^^ is s~^/^(— r)^/^ times an even function of r. 
Nonlinear scaling field analysis then predicts, in the absence of corrections, that the "scaling" 
part of the susceptibility is s~^^'^\t\'^^^ times an even function of r. Thus although our results 
for F± are not consistent with the complete absence of correction terms as discussed in section 
7 the prediction that y/sF± is even in r does appear to be preserved to all orders. 

To determine the coefficients in F± we return to the unsealed x± of sqs. (|1.9|Jl.lCl|) so 
that the terms in F± are now singularities of the function and contribute to the high order 
coefficients in the series. The cases T < and T > are again similar; for T < the 
unwanted contributions from the "short-distance" part of x are subtracted as in the s-plane 
analysis described in section |6.2| . We use essentially the same smoothing except for a 1/4 
shift in power necessitated by the difference in the singularity structure generated by the x 
rescaling of section |6.2| . That is, we replace the remainder eq. ( |6.16|) by 



n 



(6.19) 



and reduce the remainders R~ in ( |6.19D by a least squares fitting to the unknown F_ co- 
efficients in a procedure similar to that described in |TlI]. Fitting intervals An > 128 are 
typically used, and an FFT of the residuals is very useful as a diagnostic to interpret the 
observed oscillations in the residuals in terms of x singularities on the circle \s~^\ = 1. Be- 
cause the highest order terms in F_ are not fixed (i.e. known) unlike the "short-distance" 



^^See also the discussion of the quadratic recursion relations in Itzykson and Drouffc |18| where rescaling 
by was used to simplify the scaling limit. 
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terms, they tend to float and become effective amplitudes that incorporate all the higher 
order effects including the "short-distance" contributions that have not been subtracted. 
One technical result of this is that the residuals we observe in ( |6.19| ) are some six orders of 
magnitude smaller than the residuals we obtained in ( |6.16D . In part this means that whereas 
was more than adequate as a smoothing operation, ( |6.19D is marginally so. In addi- 
tion, we observe here for the first time one of the singularities from x^^\ and using the cut 
information given in eqs. (|3.26|, |3.27|) can subtract it as 



^2m ^ fs:2m + 7^15/(16rf^/2)(3V7/7r^)3/^{cos(2marccos(3/4)-7r/8) (6.20) 
- 215V7/(32m) sin(2marccos(3/4) - vr/S) + 0{l/rh?)}, 

where m = m — 1/4. 

A similar analysis has been carried out for T > and our estimates for the coefficients 
of F± are given in the Appendix. The coefficients through O(r^) are unambiguously rational 
and have been fixed in the final fittings. We have also set to zero all coefficients of odd powers 
of r in the product ^/sF± to 0(r^^) but have allowed variable coefficients of r^^, r^^ and r^^. 
These latter coefficients are, as expected, sensitive to whether we stop the "short-distance" 
subtraction at 0(r^^) or add in the additional terms we have estimated by Pade methods. 
The terms quoted in the Appendix on the other hand are completely stable and thus we 
believe reliable except possibly for the last digit. 

When we relax the constraint of zero amplitude on individual odd terms in \/sF± 
for integer 5 < /c < 15 we find no significant improvement in our fits and the resulting 
amplitudes are consistent with zero. For example, when T < T^, we find best fit coefficients 
of r'^, k odd, that in absolute magnitude are all less than 4 x 10^'^^'^^. For T > the 
corresponding bounds can be as much as 100 times larger. But in both cases these bounds are 
of a magnitude similar to what we estimate is the uncertainty in the even order coefficients 
given in the Appendix. Thus we believe \/sF± even in r to be an exact symmetry of the 
scaling-amplitude function. 

An alternative analysis using the traditional variable v = tanhi^ was also carried out. 
The natural boundary singularities at \s\ = 1 are mapped to two circles |f ± 1| = V^- 
In this expansion variable, the ferromagnetic and anti-ferromagnetic critical points are at 
V = ±(-\/2 — 1) respectively, and all other points on the two circles are farther away from 
the origin. Hence the amplitudes of any other singularities are exponentially damped and 
may be neglected in the analysis. The two analyses are in complete agreement, but the more 
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detailed s-plane analysis provides greater precision. For contemplated future analyses on 
other lattices, for which less is known about the natural boundary singularities, it may be 
necessary to use the f-plane analysis. 



7 Expected scaling form of the susceptibility. 

The basic scaling Ansatz for the singular part of the free-energy of the two-dimensional Ising 
model is 

+ gty±igh/\gt\''^'\{gu,\gtr^^''}). (7.1) 

Here gt,gh,guj are nonlinear scaling fields associated with the thermal field t, the magnetic 
field h and the irrelevant fields {uj}. The exponents yt,yh > are the thermal and magnetic 
exponents, and yj < are the irrelevant exponents. 0. The nonlinear scaling fields have 
expansions 

gt = ^a2„(t,M)/i^", ao(0,u) = 0, 

n>0 

gh = J2^2n+l{t,u)h^^^\ 

n>0 

gu, = J2c2n{t,u)h'\ (7.2) 



n>0 



where a2n, &2ri+i, C2n are smooth functions of t and u = {uj}. 

If irrelevant fields are neglected, then the known zero field free energy forces the equalities 
K|_(0) = VI (0) and ^+(0) = ^-(0). Furthermore, the absence of logarithmic terms in the 
known magnetization and the divergent part of the susceptibility requires the derivatives 
y^(0) and 1±'(0) to vanish. Aharony and Fisher have conjectured |2^, almost certainly 



correctly, that there are no logarithms multiplying the leading power law divergence of all 
higher order field derivatives, in which case the Y± are constants and analyticity on the 
critical isotherm for h ^ demands = Y_. With all these constraints built in, one can 



^^For the two-dimensional Ising model, yt = 1 and ?//,. = This scaling Ansatz assumes only a single 
resonance, between the identity and the energy. That the dimensions are integers implies higher powers of 
\ogt. The surprise is the power of t at which higher powers of logt enter. We thank Andrea Pelissetto for 
this clarification. 
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show the scahng Ansatz ( [7.1D together with the field expansions ( [7 .21) lead to 



f{t,h = 0) = -A{ao{t))Hog\aoit)\+Aoit), 
M{t<0,h = 0) = Bbi{t)\ao{t)f , (7.3) 
f3-'x±{t,h = 0) = C±(6i(t))>o(t)r^-^«2(t)ao(t)log|ao(t)| + ^(t), (7.4) 

where A, B, C± and E are constants and /5 = 1/8, 7 = 7/4. The free energy and magnetiza- 
tion eqs. ([TT^ ) determine the scaling field coefficients ao(t) and bi(t) in (|7.2|) . The presence 
of irrelevant scaling fields will be expected to manifest themselves as deviations from the 
predicted form of the susceptibility in (7.4) and/or as deviations from the unique prediction 
for the coefficient of C±. 

Working as usual in the temperature variable r = (1/s — s)/2, we write the predicted 
isotropic susceptibility from (7.4) as 

p-'x±{r,h = 0) = Coj2ir,y2)7/4|r|-^/4F(A&F) -i?o/(2i^cv^)rlog|r|eo(r) + Do(r) 

(7.5) 



where F(A&F) has already been given in eq. ( |1.19 ). The "short-distance" contribution to 



X is here predicted to be given by the sum of the term containing 69(1"), which arises as the 
mixing of the first two terms in the expansion of gt in ([7.2|), and the analytic Dq{t). 



The clear implication of eq. ( 1.18 ), which shows the exact F± is not equal to F(A&F) 
is that irrelevant variables do play a role, and the multipliers F± represent the contribution 
of a number of scaling fields. While there are suggestions in the literature for what these 
scaling fields might beQ it is unlikely that a unique identification could be made here since 
we are dealing with a single isolated model with no free parameters to vary. A corresponding 
analysis of the anisotropic square, and the triangular and hexagonal lattices is likely to be 



enlightening in this regard. In |^ a study of difference equations is given, which implicitly 



outlines what is needed to obtain the difference equations for the hexagonal and triangular 



lattices. Additional material in this respect can also be found in ||50|] 



An excellent early discussion of the types of corrections that might be found, together with a search for 
some of them can be found in Blote and den Nijs Developments in our understanding of the predictions 
of conformal field theory |^-^ lead us to believe that a virtually complete explanation of corrections to 
scaling is obtainable, at least in principle. Our analysis supports the conclusion of Barma and Fisher that 
a correction-to-scaling term with exponent 61 = 4/3 is absent for the pure S ~ 1/2 Ising model susceptibility 
considered here. A mechanism in terms of generators of the energy family of the Virasoro algebra is adduced 
by Caselle et al. Q which gives rise to corrections at order t'^ as observed. 
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Note however that two different effects manifest themselves. At fourth order in t (or, 
equivalently, r) scahng under the assumption of only two nonlinear scaling fields breaks down, 
as evidenced by the difference between the coefficients of in eqs. ( |1.18| ) and ( |1.19| ). However 
the corresponding high- and low-temperature amplitudes still satisfy Cq /Cq = /C~ for 
all j < 5. For j > 5, not only does simple scaling fail to hold, but this equality also breaks 
down. 

In the vicinity of the anti-ferromagnetic point in the high-temperature phase, x is given 
exclusively by a "short-distance" term (|1.16| ). Both i?af and Bf from eqs. ( |1.16D and ( |1.17D 
have expansions of the same form (|1.15|) , where the sum over p is restricted to < q. The 
coefficients in this expansion can be determined from the short-distance correlations, and 
accurate values for the expansions of -Bf/af are given in the Appendix through 0(r^^). 

Again we note that there are terms in these "short-distance" functions that are not of the 
Aharony and Fisher [0,|2^ predicted form (|7.5|) based on the absence of irrelevant variables. 

We conclude with some speculative remarks. We cannot account physically for terms of 
order t''(log with p >2 and q > p"^ in Bf, though we can see their origin mathematically, 
as discussed below eq. ( [4. IS] ) and in footnote 12. While terms without logarithms and terms 
of order ^(log |t|) are expected, it is surprising (to us) that higher powers of log \ t\ enter at 
the orders they do. 

From conformal field theory, we have predictions for the irrelevant exponents yj = 
—2, —4, —6, .... The fact that f^^ = f^^ for k < Q, though f^^ is not equal to the corre- 
sponding term in (|1.19| ) for k = 4,5 suggests the presence of only a single irrelevant operator 
contributing at order r^, while the breakdown of high-low temperature symmetry in F± at 
O(r^) suggests that more than one scaling operator couples to the lattice magnetization at 
this order. A corresponding study to that reported here on the triangular and honeycomb 
lattices, as well as on the anisotropic square lattice is likely to be enlightening, and we hope 
to report on this in future. 
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Appendix 



Here we list the expansion coefficients bp^q of the "short-distance" functions -Bf/af defined 
in ( |1.15| ) to 0(r^^). The prefactor (vT+T^+r)^/^ = is understood to be expanded in a 

series in r. The leading constant and coefficient of r log |r| have been reported previously 
The F± series as deduced in section ^]3] are also given. 



[-0.104133245093831026452160126860473433716236727314 

- 0.07436886975320708001995859169799500328047632028r 

- 0.0081447139091195995371542858655723893266057740r2 
+ 0.004504107712232015926355020852986970591364528x3 
+ 0.23961879425472180967837072450742931180586109r^ 

- 0.0025399505953392329612162686121616238176205r^ 

- 0.235288909669962491804066210882350821445764r'^ 
+ 0.00191570753170091141409998516033460855797r^ 
+ 0.2143400966115384518711435705343125612378r^ 

- 0.000883215706003328768611915486246075323r^ 

- 0.19422062840719623752953468278129284679r^° 

+ 0.0000072335097772632765778839359680528r^^ 
+ 0.177102037555467190714704023746648559r^2 

+ 0.00068881 10962684387331860926084517r^3 

- 0.16279253648974618861881216566686r^^ 
+ log |r| 

(0.032352268477309406090656526721221666637730948898r 

- 0.0057755293796884630091487564013201013677152980r3 
+ 0.059074961290345476578516085774495545264759330r^ 
+ 0.00305849157585622544005057759535229287938174r^ 

- 0.0591662722088409053375931018028970139567911r^ 
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- 0.002067088393167114141650194740281136875636r' 
+ 0.05424693070421409615112542698595864778919t^ 
+ 0.0010601025315498815900774416057541651837t'^ 

- 0.049300253157082567741316861339709144063t^° 

- 0.00026830064161204203467706137637400358t^^ 
+ 0.0450270525719569186212816103126356308r^^ 

- 0.00034332683257234543036792535081332r^^ 

- 0.041428586463052869356803144137620T^^) 
+ (log|r|f 

(0.0093915698711458721317953318727075770649513654T^ 

- 0.00869592546287923802156416645191752987912922T^ 
+ 0.007669481493104540876445085447422616885330r^ 

+ 0.00015428438297902275440225213783285077606r^ 

- 0.0068054076881441249098452112921129773269T^° 

- 0.000310520937481414524040686012525223279T^^ 
+ 0.00611386643219454473116391019937140965x^2 
+ 0.000444606198235804033861443998682830t^^ 

- 0.0055571002151161308034896964314679r^^) 
+ (log|r|)3 

(-0.000015771569138451840480001012621461738178r^ 
+ 0.0000344282066208887553647799856857753380t^^ 

- 0.0000524427177487226174161583779149393t^^)] , 

[0.15886652296094748823335923X3690210116925239008416 
+ 0.149566836938535905194382029433591286374711207262T 
+ 0.01071222587983288033470968550659996768542030678t2 
+ 0.012753018839962401953955207805215360913467497lT^ 
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- 0.01 1741 188869656263932121387296300743594029390r^ 
-0.01406604087566590060620992322775625815515533t^ 
+ 0.0131064546156258402249424759665220798848681t^ 
+ 0.012239696625538370626786005459530159711716t^ 
-0.01184019404541084813958321002995560636877t^ 
-0.0105854093023116661362507232392645147231r^ 

+ 0.010151560037724359473553197335905170854r^° 
+ 0.0090800041123311937161054945314071828lT^^ 
-0.0085420122287896456879459087815054030r^^ 
-0.00771702694013238358077176900242074t^^ 
+ 0.007123677682511208149032476379667t^^ 
+ log |r| 

(-0.1553171901580110585934133538932734529992121600305r 
+ 0.03206714814586975221843437287457551882247161782r3 
-0.0077168875724615093064542922101962689299768599t^ 
-0.015675211573817078943430169665269657287132910t^ 
-0.00028554245153720354627897919087710530890677t^ 
+ 0.0096072545027321808179041903535130201897217r^ 
+ 0.004835406420625092236673413378307375358908r^ 
-0.00606499034448050751379194815071149626812t^ 

- 0.007340015041447402356245461 1060746360875t^° 
+ 0.003910356521403913091050321141297009252t^^ 
+ 0.00870842744568158003036434762719697635t^^ 
-0.002697783010884752101384006121375890t^^ 
-0.0094056230380765607719474925088649r^^) 

+ (log|r|)2 

(0.01153371437882328027949011442761203640684043805t^ 
-0.011311734920691560067535056532207842716405684t^ 
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+ 0.0100457687111988577404299867962466051265974r« 
-0.000475698571097159420906450182271928179428r^ 
- 0.00878397202228689639470985683437717938463r^° 
+ 0.0011571801729636538264100914359800686355r^^ 
+ 0.007680651 109512704070606639646988801296r^2 

-0.0018650912616201532939412831153215046r^3 
-0.00674470189451526288478200059343432r^^) 
+ (}og\r\f 

(0.0000578997194764877297760067221144062249541r^ 
-0.00016991508824012890240796446744935908812r^^ 
+ 0.00032664884687465587957270016883093909r^3)] 

F+ = 1 + r/2 + 5rV8 + 3rVl6 - 23rV384 - 35rV768 - 0.1329693327418753330r^ 

- 0.05899768720427100r^ + 0.121586869804903r^ + 0.0766007994119r^ 

- 0.10751871874r^° - 0.078346589r" + 0.0960583r^2 ^ o.07592r^3 

- 0.087r^^ - O.lr^^ + . . . 

= (Vl + r2 + r)^/2(l + - rVl2 - 0.1235292285752086663r'^ 

+ 0.136610949809095r^ - 0.13043897213r^° + 0.1215129r^2 _ o.ll3r^^ + . . . ) 

F_ = 1 + r/2 + 5rV8 + 3x^16 - 23rV384 - 35rV768 - 6.330746944662603289734r^ 

- 3.1578864931646349782r^ + 5.462251 18896595954r*^ + 3.521655160482472r^ 

- 4.6602157191837r^° - 3.40963923001r" + 4.055875878x^2 ^ 3.2008085r^3 

- 3.59746r^^ - 2.985r^^ + . . . 

= (Vl + r2 + r)^^^(l + rV2 - r^U - 6.321306840495936623067x6 

+ 6.25199747046024329r^ - 5.6896599756180r^° + 5.142218271r^2 

- 4.67472r^^ + . . . ) 

The last digit in each term above may not be rehable. 

For completeness we give also the leading susceptibility amplitudes evaluated to higher ac- 
curacy than reported in [Q: 

Co+ = 1.000815260440212647119476363047210236937534925597789(2is:eV^)-^/^V2 

Cq = 1.000960328725262189480934955172097320572505951770117(2i\:cV2)"^/^V2/(127r). 
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